Low-energy sub-gap states in multi-channel Majorana wires 
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One-dimensional p-wave superconductors are known to harbor Majorana bound states at their ends. Super- 
conducting wires with a finite width W may have fermionic subgap states in addition to possible Majorana end 
states. While they do not necessarily inhibit the use of Majorana end states for topological computation, these 
subgap states can obscure the identification of a topological phase through a density-of-states measurement. 
We present two simple models to describe low-energy fermionic subgap states. If the wire's width W is much 
smaller than the superconductor coherence length ^, the relevant subgap states are localized near the ends of the 
wire and cluster near zero energy, whereas the lowest-energy subgap states are delocalized if TV > f . Notably, 
the energy of the lowest-lying fermionic subgap state (if present at all) has a maximum for ~ ^. 
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The search for Majorana fermions has attracted a great 
deal of interest in the last few years |1|. Notably their non- 
local properties and non-abelian braiding statistics make Ma- 
jorana fermion systems attractive candidates for fault tolerant 
quantum computation |2-4|. The present wave of interest is 
driven by a number of proposals that suggest ways of realizing 
and manipulating Majorana states in solid state systems, most 
prominently interfaces of s-wave superconductors and topo- 
logical insulators fSjISJ, half-metallic ferromagnets lITHll, or 
semiconductor films or wires liTOl - frSl . where the latter stand 
out because Majorana manipulation require a mere series of 
gate operations ifTSll . In all these proposals, the proximity cou- 
pling to the s-wave superconductor effectively turns the nor- 
mal metal into a p-wave superconductor, which is well known 
to harbor Majorana fermions at its ends or edges lfT4Hr7l . 

Majorana bound states at ends of what is effectively a p- 
wave superconducting wire can be analyzed most straightfor- 
wardly if these wires are strictly one dimensional, with only a 
single propagating mode at the Fermi level in the absence of 
superconductivity iTO","!!!- Nevertheless, Majorana end states 
can also exist in a quasi one-dimensional geometry. The effect 
of multiple transverse channels, present in most realistic real- 
izations, has been addressed in Refs. lfT2l [1814211 . Specifically 
one sees that a complex p-wave superconductor in a strip ge- 
ometry undergoes a series of oscillatory quantum phase tran- 
sitions between topologically trivial and topologically non- 
trivial phases (without and with Majorana end states, respec- 
tively) as the strip width W or chemical potential are varied. 
Both with and without Majorana end states, a range of sub- 
gap states is found [19J, analogous to the sub-gap states in 
vortex cores of bulk superconductors ll22l . Although the mere 
presence of sub-gap states does not prohibit the use of Majo- 
rana end states for topological quantum computation the 
presence of low-lying sub-gap states clearly obstructs an un- 
ambiguous experimental verification of the Majorana states. 

The purpose of this paper is to systematically analyse the 
energies of the sub-gap states as a function of the sample 
width W . Throughout we assume that ^ > Ap, where 
Ap = 2Ti/k-p is the Fermi wavelength in the absence of su- 
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FIG. 1. (Color online) Excitation spectrum of a p-wave supercon- 
ducting strip as a function of the strip width W , for small Wkp. The 
solid (red) curves and the shaded region indicate the discrete energies 
of the end states and the continuous spectrum, respectively. The data 
shown in the figure is obtained from a solution of the effective mod- 
els described in the text with fcp^ = 100. Data points obtained from 
a numerical solution of the lattice with — 4t < fi < —3t are indistin- 
guishable from the curves shown here. Inset: Density of fermionic 
end states in the limit Ap < W < ^. Here eiyp = tvA{W/^)^ and 

l^typ = Wkp/iTVetyp). 



perconductivity and ^ the superconducting coherence length. 
Our results, which will be discussed in detail below, reveal 
two different regimes, depending on the relative magnitude of 
W and 

I. lfW<£, there is an alternation of topological and non- 
topological phases, where each phase harbors at most 
one Majorana state at zero energy and Nf < N ^ 
int [iW/T:)y^kp — ^^^J fermionic subgap states local- 
ized near each end of the wire, see Fig. [T] The lowest- 
lying fermionic subgap state (if present) has energy 



^min '—^ ^2 



WXt, 



^2 \n{W/\f 



(1) 



one has N{ = 



For wire widths W ^ 
N and the highest-lying subgap state has energy 
AW^/£,'^Xp, well below the bulk excitation gap A. 
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2. For W > there is an alternation of topological and 
non-topological phases, with Majorana end states in the 
topological phase, and delocalized subgap states with 
minimum energy 



Below we will set 



A,, and formulate our main re- 



e • -Ae-^/« 
Equation (j2|i was obtained previously in Ref.lZTl 



(2) 



We note that as most experimental proposals involve a rather 
weak induced superconductivity because of the presence of 
Schottky barrier or the dependence on spin-orbit coupling, we 
believe that the quasi one-dimensional regime W/^ <C 1 is 
most relevant for possible applications. Hence, these find- 
ing imply that, in the experimentally accessible regime, local- 
ized subgap states cluster around zero energy, thus obstruct- 
ing the experimental verification of Majorana states through a 
density-of-states measurement, see Fig.[T] In fact, the optimal 
wire width W that gives rise to the largest excitation gap to 
fermionic subgap states is ^ ^ (Fig. |2]i, although the op- 
timal excitation gap for a multi-channel wire is always lower 
than that of a single channel wire. Nonetheless, for W/ C ^ 1 
there is no coupling between subgap states localized at op- 
posite ends of the strip, as the excitation gap to delocalized 
bulk states is large (> A, see Fig. [T]i, making the narrow 
limit favorable for topological protection. Our analysis pro- 
vides a straightforward physical interpretation of the oscilla- 
tory quantum phase transitions and highlights broad similari- 
ties between the Z to Z2 crossover and the nucleation of topo- 
logical phases discussed for example in Il24ll25l . 

For the detailed discussion of these findings we employ a 
two-dimensional lattice model of a complex p-wave supercon- 
ductor ifTSl [TTl [26i . The same model was studied recently in 
Refs. |[T8l[l9l|27l|28l. We consider a strip geometry of length 
L = N^a and width W = Nya, where a is the lattice con- 
stant, with the Hamiltonian 

1=1 j=l 



E 

i=i j=i 



h.c 



J = l 

E E (4,+iCM +A,4.+i4.+h.c.), (3) 

i=i j=i 

where fi is the chemical potential, t the hopping amplitude, 
and Aj; and Aj, set the p-wave pairing amplitudes in the lon- 
gitudinal and transverse directions, respectively. We will fo- 
cus on the scenario that the p-wave superconducting order is 
inherited from proximity coupling to a bulk superconductor, 
as in the proposals of Refs.^HISl so that no self-consistency 
conditions need to be accounted for The bulk dispersion of 
Eq. Q gives excitations at energies 



El 



[—fi — 2t(cos fcjjfl + cos kya)]'^ 
+ AAl sin^ kj^a + 4A^ sin^ kyU. 



suits in the continuum limit, which is obtained by send- 
ing a while keeping L, W, and the number of elec- 
trons fixed. In this limit, /i — > —4t + t{kpa)^ and — > 



where fcp = 

,2 



27r/A 



F IS 



kpWC^^ ~ kDHfivp)^ + fc2A2, 

the Fermi wavenumber, hvp — 2tkpa'' the Fermi veloc- 
ity, and A — 2Axkpa — 2Aykpa the bulk excitation gap. 
The superconductor coherence length is ^ = hvp / A. The 
model ([3]) falls into class D in the symmetry classification 
of non-interacting topological insulators and superconductors 
II29I l30i . In one dimension, the archetypal model in the D 
class is the so called Majorana chain or wire ifTSl . which is 
characterized by a Z2 topological invariant. In the two dimen- 
sional (2d) limit, class D is characterized by a Z topological 
invariant. This model, and its — iky partner correspond to 
realizations with Chern numbers of cr = 0,±1; examples of 
higher Chern numbers do exist, see e.g., | fl6ll3T1433l . 

For W < 7r(fcp — ,^2)-i/2 (jjg model Oil is a trivial insula- 
tor Upon increasing the sample width, the system undergoes 
a series of topological phase transitions. In order to study low- 
lying excitations in each phase, we use two separate calcula- 
tions, valid for <C ^ and W > respectively. 

The case <C ^. If ^ ^, hypothetical end states have a 
vanishing expectation value of the transverse momentum py. 
For this reason, it is a good approximation to treat the term 
proportional to A^ in the Hamitonian ( 31 in perturbation the- 
ory. Solving for the eigenstates of Eq. ( 3JI with Aj, = in the 
limit of large L, we find (in the continuum limit) 



N = int [{W/tt) 



kp 



(5) 



zero-energy Majorana states at the left (+) or right (— ) ends 
of the wire, with creation/annihilation operators jn.± — 
e^'"I^Y.i j Un,±iia,ja){ctj ± icl n 1, 2, . . . , iV, and 



Un,+ ix,y) 



Here k„ 



kp 



n-^/^ sin(n7ry/W^) sin(A:„^a;)e~="/«. 

(6) 

(mr/W)'^ — fin is a normalization 
constant, and 7i„._(a;, y) = u„ a;, y). Full results for the 
lattice model are qualitatively similar f34\. The pairing term 
proportional to Ay is then treated in degenerate perturbation 
theory. Using the N zero-energy states of Eq. (|6]l at each end 
as a basis, one finds the effective N x N Hamiltonian Hnm 
with Hnm = if n + m is even and 



32iAXpmn 



(7) 



A|n2)(W2 



A2m2) 



8A2(n2+TO2) + M^(„2 



to2)^ 



(4) 



if n + m is odd. The eigenvalues e of the antisymmatric ma- 
trix (|7]i come in pairs ±e; the Hamiltonian (|7]l has a single 
zero eigenvalue if TV is odd. In that case the corresponding 
zero mode at the left/right end is a linear combination of or- 
bitals u„ -t (x, y) with n odd only. The topologically trivial 
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phases, which have no zero-energy end states, occur for even 
N. Since N increases stepwise as the chemical potential /i 
or the width W is varied, the system thus undergoes a se- 
quence of topological phase transitions. We could not obtain 
a closed-form expression for the eigenvalues of the effective 
Hamiltonian (|7| for general N. A simple scaling analysis 
of the Hamiltonian (j7]l shows that the lowest positive eigen- 
value e„iin is of order AAp min(W^/^^, 1/W), so that the op- 
timal (largest) separation between the Majorana state and the 
lowest-lying fermionic excitation is achieved for W of order 
^. A numerical evaluation of the spectrum reveals an addi- 
tional logarithmic correction in the limit W <^ thus giving 
the estimate ([T]i for ^ ^. The highest eigenvalue of the 
matrix (j?} is of order Amax(l, VF^/Ap^^). For < f , the 
density of states has median Styp ^ 7rA(VF/^)^, well below 
the bulk excitation gap A, thus justifying the perturbative pro- 
cedure, see Fig. [T] (inset). 

For W ^ ^, the spectrum of end states found here is similar 
to that of the subgap states in a vortex core ll22l . which appear 
at a regular spacing ^ AAp/^, filling up the entire region be- 
tween zero excitation energy and the bulk gap A. On the other 
hand, in the limit that the wire width is small in comparison 
to the coherence length, the end states cluster near zero exci- 
tation energy and their spacing is anomalously small. 

The case W > The end states discussed sofar are not the 
only possible low-energy excitations of the system. In addi- 
tion to the end states, there are lateral Majorana modes local- 
ized near the edges at y = and y — W llT4l[T6ll . At transition 
points between the topological phases, the lowest energies of 
these lateral edge modes drops to zero — see the discussion 
below Eq. ([8]) for details — , so that the effective Hamilto- 
nian (|7]i provides a complete description of the low-energy 
excitations away from the transition points only. Moreover, 
away from the transition points, the energy splitting of the 
lateral modes decreases exponentially as the width increases, 
and eventually they become the dominant low-energy excita- 
tion for sufficiently large W > with an excitation gap that 
follows - Ae-^/«. 

In order to gain a deeper understanding of the role of the lat- 
eral Majorana modes, we now apply periodic boundary condi- 
tions in the longitudinal (x) direction, while keeping the hard 
boundaries at y = (+) and y — W (— ). The calculation 
presented here allows one to calculate the excitation gap to 
the lateral edge modes in the continuum limit and reveals that 
the Z2 topological order of the bulk can be understood as ef- 
fective Majorana hopping model. The Majorana modes on the 
lateral edges are labeled by the wavenumber and have op- 
erators 7fe^,± = E™,„Wfe.,±("a)e''''"""'(cm,„ ± c3„,„) and 
energy e±{kx)- In the limit W/£, — > 00 the wavefunctions 
Uk^,± {y) are known lfT6l . 



(2/) 



sm{kyy) 



ky — 



2,2 _ 7,2 _ 



(8) 



iuk^^+{W — y), and their energy is 

e±{kx)^±akx), ^ikx) = kxA/kF. (9) 

The wavefunctions ([8]l also satisfy the boundary conditions 
for finite W if kyW/ir is an integer These special points, 
at which the spectrum of the lateral Majorana modes is gap- 
less, correspond to the boundaries between the topological 
and non-topological phases |21|. [Note the consistency with 
Eq. (|5]l above.] In the vicinity of the phase boundary at 
ky = ky — niyTT/W, niy integer, the energies of the lat- 
eral Majorana modes can be obtained in perturbation theory, 
which leads to the effective 2x2 matrix Hamiltonian 



^(kx) h{kx) 

HK) -i{kx) 



for the lateral Majorana modes, with 



Kkx)^-^^WyY 



kt 



(10) 



(11) 



with rtk^ a normalization constant and Uk^-{y) 



where the superscript "c" refers to the wavefunctions at the 
transition point, i.e., with ky — ky. The eigenvalues of H^^ 
are e±{kx) = ±[/i(fca;)^ + Hkx)'^]^^'^ ■ Interestingly, the sign 
of h{kx) encodes the sign of the effective mass of the sys- 
tem lfT4ll and therefore also encodes the Z2 topological in- 
dex of the system ifTSl [T6l . Although the effective Hamilto- 
nian ( lOi has been derived in the vicinity of transition points 
between topological and non-topological phases {\ky — ky] 
small), comparison with a full numerical solution of the lat- 
tice model Q shows good agreement for all W if my is taken 
to be the integer closest to kyW/n. (Agreement within 13% 
for ^ ^, data not shown.) Making use of the estimate 
Kx.+ l^tfcx,-> ~ {-lY'+'^{W/2Cje-^/^ for W > 
find that the typical gap e„iin of the lateral Majorana modes 
away from the transition points is given by Eq. (|2]i if ^ ^. 

Figure |2] shows the value of the energy e„iin of the lowest- 
lying fermionic excitation at the mid-point of each topological 
phase, as a function of the width W . We observe that the two 
calculations outlined above are in excellent agreement with 
the full numerical calculation in the limits W ^ and W ^ 
^. Although there are quantitative deviations in the crossover 
between these two limits, £„ii,i is qualitatively well described 
by the minimum of the two asymptotic expressions (|2]) and 

Conclusion and perspective - We have presented two 
simple analytical models to describe two types of low-lying 
fermionic excitations in a multichannel p wave superconduct- 
ing wire. The first model describes fermionic end states, 
which are the lowest relevant excitations for wire widths W 
well below the superconductor coherence length ^. The sec- 
ond model applies to delocalized excitations along the wire's 
edges, which are the relevant low-energy excitations if > 
^. While our results were obtained for a simplified model of a 
pH-ip superconductor, we expect that the general scaling prop- 
erties presented here will remain present in a more realistic 
treatment of spin-orbit coupled superconductors. 
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FIG. 2. (Color online) The energy Smin of the lowest-lying fermionic 
excitation in the topological non trivial phases, as a function of the 
sample width W . The main panel shows the value of Emin at the mid- 
point of each topological non-trivial phase for /i = — 3t (see inset), 
which illustrates the large-scale dependence of e-mi-n on the sample 
width W . Squares and triangles refer to full numerical solutions of 
the lattice model ijsjl for kpS^ — 50 and 33| respectively; The dotted 
solid black curves and the red solid curve represent solutions of the 
effective models (jTJ and \\Q) , respectively. Inset: Emm vs. W for the 
topologically non-trivial phases for /cf^ = 50, obtained by combin- 
ing data from a numerical solution of the lattice model |3j for various 
values of /i near —St. 



The fermionic end states, which are the relevant low-energy 
excitations for VF ^ ^, appear in the topological as well as 
in the non-topological regimes, as soon as the wire has more 
than one channel, W > TT{kp — ^^^)^^/^. For thin wires, 
<C ^, the fermionic end states cluster around zero energy, 
see Fig. [T] so that the mere observation of a large local den- 
sity of states near zero energy at the end of the wire is not a 
reliable signature of a Majorana bound state. Energy differ- 
ences of order AApVF/^'^ must be resolved if one wants to 
separate the fermionic end states and the Majorana end states. 
This is a stricter requirement than for energy levels in vortex 
cores, which are spaced by a distance ^ A{Xf/£,) ll22l . The 
distance Emin to the lowest-lying fermionic subgap excitation 
has a maximum for wire width W ^ The optimal value, 
£min ~ AAp/C^ has the same scaling as in a vortex core. 

While the end states obscure the identification of a topolog- 
ical phase through a density-of-states measurement, the pres- 
ence of localized fermionic states near the wire ends does not 
necessarily inhibit topological quantum computation based on 
the Majorana states ll23l |35]| . In fact, the large gap ^ A 
to the delocalized bulk excitations, makes the narrow limit 
W/£, ^ 1 favorable for topological protection. On the other 
hand, the lateral Majorana modes are delocalized, and their 
presence at low energy for W > poses a restriction on the 
use of the Majorana end states for topological quantum com- 
putation. 

It is important to point out that the end states we discussed 
here exist at the ends of multichannel superconducting wires. 
The situation is different at an interface between two regions 
in which the channel numbers changes by a small amount 
only. This is the relevant scenario if for semiconductor-based 
topological superconductors, where the change in channel 
number is induced by a small change in a gate voltage. In 



that case the number of fermionic subgap states localized at 
the interface is at most half the difference in channel number 
In particular, if the channel number changes by one, the inter- 
face harbors a Majorana fermion, but no low-lying fermionic 
states. 

We close with a remark on the effect of disorder. Even weak 
disorder is known to lead to the formation of subgap states lo- 
calized somewhere in the bulk of the superconductor, with an 
algebraic density of states low energy 1 17, 36 1. The effect of 
weak disorder on the low-energy end states, in contrast, is ex- 
pected to be small, because the disorder potential has vanish- 
ing matrix elements between the basis states (j6]) that are used 
for the construction of the effective end-state Hamiltonian (|7]i. 
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Appendix A 

The bulk dispersion relation in two dimensions is given as 
Ek - v/^^ + lap with 

hk = —jJ- — 2t(cos kxO. + cos kyo), (12) 
a = 2A(sin fej-a + i sin fcj,a). (13) 

The ground state may be described by the BCS state lfT6l 

|gs)= {uk + Vka^y,a}_y,)\v2ic) . (14) 

(k-k) 



where the product is over distinct (fc, — fc) pairs and 

Uk = ^Jl/2{l + hklEk), 

Vl/2(l-/^fe/i^fe)e"^^^^ 



1,2,. 



. , N, and 



(15) 
(16) 



This model has a topological phase transition at — At. The 
phase with > At corresponds to the strong pairing phase 
of Read and Green lfT6l and the Toric Code phase of the vor- 
tex free Kitaev honeycomb system [37 1. When < At the 
model is in the non-abelian weak pairing phase of Read and 
Green that supports bound Majorana fermions at vortex cores 
and along domain walls separating regions of topologically 
trivial and non-trivial phase. This phase is closely related with 
the Pfaffian filling fraction = 5/2 fractional quantum hall 
state and the gapped B-phase of the Kitaev honeycomb system 

lElllTl. 



Appendix B 

We summarize the main results of the finite-width strip in 
the limit W ^ for the lattice model. Solving for the eigen- 
states of Eq. (3) in the main text for Aj^ = we find 



N = int 



1 



arccos 



-2^ 



2t 



zero energy IVIajorana states with creation/annihilation oper- 



ators Jn,± 



J2i.Un,±{ia,ja)ic,^j ± ict -), 



Un,+ (x,y) 



Here 



Ny + 1 " 



f^-VV/'^sinf^^ 



(18) 



(19) 



Un,-{x, y) = Un.+ {L—x, y), and is given by an analogous 
expression to ftk^ see Eq. ([32]). Accounting for the pairing Ay 
in degenerate perturbation theory then results in the effective 
N X N matrix Hamiltonian 




Hri.vn. — 



AiAy sin fc,„a sin fc„a 
{Ny + 1) COS femO — cos kna 

4AV(^l-4t2)(62 _4t2) 



where 6„ = fj, + 2t cos fc„a. In the continuum limit, one re- 
covers Eq. (7) of the main text. 



Appendix C 

Here we review the calculation of the edge mode structure 
for the two-dimensional strip with periodic boundary condi- 
tions in the longitudinal (x) direction. If the boundaries at 
y — (+) and y ~ W (— ) are far apart, it is possible to 
exactly solve for the edge spectrum and recursively resolve 
the corresponding edge modes for arbitrary values of k^- It 
is well known, see for example ||T6l |38l that the modes on 
the well separated (upper or lower) edges are described by the 
Majorana ansatz 



= E "fc.,±('^«)e''^'"''(cL,„ ± c„,„) (21) 



For the lattice model in question, the bottom (top) edge spec- 
trum is given by i? = ±2| Aj; | sin(fcj;a), see for example Refs. 
ll2n l39l l. Using this energy eigenvalue we can, on the bot- 
tom edge for example, calculate the unnormalized values of 
Uk^.+ ina) recursively through the transfer matrix as 



with 



"fcx,+ ((" + l)a) 



n 



Uk^^+{{n - l)a) 
Uk^^+{na) 



1 

'C/A ~B/A 



(22) 



(23) 



(17) and A = t + \Ay\, B = fi + 2tcos{kxa), SLud C ^ t ~ \Ay 



More generally we have 



Uk^,+ i{n + l)a) 



= \Tk. 



Uk^,+ ia) 



(24) 
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The eigenvalues of the transfer matrix are A± = (— B ± 
\/4AC-B2)/2A The corresponding eigenvectors are 



V± 



1 1 



(25) 



where d = \/\ + C /A. Note that these eigenvectors are not 
orthogonal. The inverse of the eigenvector matrix is 



VZ^ = — 



A_ -1 
-A+ 1 



(26) 



which allows one to write 



XI 
A" 



(27) 



1 



A_ - A. 



A_A 



n+1 



A"+'A, A 



Xl_ - xi 

n+l 



For hard wall boundary conditions where Uk^.+ {0) = we 
would have 



rpn 







1 



A_ - A, 



X"l - XI 



and by setting Ufe^, ,+ (a) = 1 we find 



XI -XI , , W^sinfcy?; 



Here 



A_ - A_ 



VcjA 



(28) 

(29) 
(30) 



1 I + 2t cos k^a) 
— arccos < ^ 



Z"2 _ 1,2 _ f-2 



(31) 



The normalization of the wavefunction is given by -v/ilfc^ 
where ilfc^ 



S2 = 



S3 



(si - S2 - S3)/2 with 

= r2 - 1 
exp( 2i{Ny + l)fcya)r2(^«+i) - 1 

(r^) exp(2i/cj,a) — 1 
exp{-2i{Ny + l)fcya)r2(^«+i) - 1 
(r^) exp(— 2ifcj,a) — 1 



(32) 



The coupling between Majorana modes at y = (+) and 
y — W given by an effective 1-d, low-energy Hamilto- 

nian, which in the Majorana basis 7fe^,± takes the form: 



where 



^fcx 



£,k^ — 2Asinfcj;a > k^A/kp 
hk, = 5e^{ul Ay)\ul iy)) 



(33) 



(34) 
(35) 



with 



a->-0 



/i + 2t cos(fca;a) + 2-\/ — |Ap cos fcy( 

.hv-p 



(36) 



and 



(37) 



(_l)™+i^e-^/« (38) 
2^ 



The energy of ( 33 1 take the usual form Ek^ 
and vanishes foffta; = when 



ky = arccos ■ 



-(/i + 2t) 



(39) 
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